home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / libcruft / blas / dsyr2.f < prev    next >
Text File  |  1997-06-25  |  8KB  |  234 lines
  1. *
  2. ************************************************************************
  3. *
  4.       SUBROUTINE DSYR2 ( UPLO, N, ALPHA, X, INCX, Y, INCY, A, LDA )
  5. *     .. Scalar Arguments ..
  6.       DOUBLE PRECISION   ALPHA
  7.       INTEGER            INCX, INCY, LDA, N
  8.       CHARACTER*1        UPLO
  9. *     .. Array Arguments ..
  10.       DOUBLE PRECISION   A( LDA, * ), X( * ), Y( * )
  11. *     ..
  12. *
  13. *  Purpose
  14. *  =======
  15. *
  16. *  DSYR2  performs the symmetric rank 2 operation
  17. *
  18. *     A := alpha*x*y' + alpha*y*x' + A,
  19. *
  20. *  where alpha is a scalar, x and y are n element vectors and A is an n
  21. *  by n symmetric matrix.
  22. *
  23. *  Parameters
  24. *  ==========
  25. *
  26. *  UPLO   - CHARACTER*1.
  27. *           On entry, UPLO specifies whether the upper or lower
  28. *           triangular part of the array A is to be referenced as
  29. *           follows:
  30. *
  31. *              UPLO = 'U' or 'u'   Only the upper triangular part of A
  32. *                                  is to be referenced.
  33. *
  34. *              UPLO = 'L' or 'l'   Only the lower triangular part of A
  35. *                                  is to be referenced.
  36. *
  37. *           Unchanged on exit.
  38. *
  39. *  N      - INTEGER.
  40. *           On entry, N specifies the order of the matrix A.
  41. *           N must be at least zero.
  42. *           Unchanged on exit.
  43. *
  44. *  ALPHA  - DOUBLE PRECISION.
  45. *           On entry, ALPHA specifies the scalar alpha.
  46. *           Unchanged on exit.
  47. *
  48. *  X      - DOUBLE PRECISION array of dimension at least
  49. *           ( 1 + ( n - 1 )*abs( INCX ) ).
  50. *           Before entry, the incremented array X must contain the n
  51. *           element vector x.
  52. *           Unchanged on exit.
  53. *
  54. *  INCX   - INTEGER.
  55. *           On entry, INCX specifies the increment for the elements of
  56. *           X. INCX must not be zero.
  57. *           Unchanged on exit.
  58. *
  59. *  Y      - DOUBLE PRECISION array of dimension at least
  60. *           ( 1 + ( n - 1 )*abs( INCY ) ).
  61. *           Before entry, the incremented array Y must contain the n
  62. *           element vector y.
  63. *           Unchanged on exit.
  64. *
  65. *  INCY   - INTEGER.
  66. *           On entry, INCY specifies the increment for the elements of
  67. *           Y. INCY must not be zero.
  68. *           Unchanged on exit.
  69. *
  70. *  A      - DOUBLE PRECISION array of DIMENSION ( LDA, n ).
  71. *           Before entry with  UPLO = 'U' or 'u', the leading n by n
  72. *           upper triangular part of the array A must contain the upper
  73. *           triangular part of the symmetric matrix and the strictly
  74. *           lower triangular part of A is not referenced. On exit, the
  75. *           upper triangular part of the array A is overwritten by the
  76. *           upper triangular part of the updated matrix.
  77. *           Before entry with UPLO = 'L' or 'l', the leading n by n
  78. *           lower triangular part of the array A must contain the lower
  79. *           triangular part of the symmetric matrix and the strictly
  80. *           upper triangular part of A is not referenced. On exit, the
  81. *           lower triangular part of the array A is overwritten by the
  82. *           lower triangular part of the updated matrix.
  83. *
  84. *  LDA    - INTEGER.
  85. *           On entry, LDA specifies the first dimension of A as declared
  86. *           in the calling (sub) program. LDA must be at least
  87. *           max( 1, n ).
  88. *           Unchanged on exit.
  89. *
  90. *
  91. *  Level 2 Blas routine.
  92. *
  93. *  -- Written on 22-October-1986.
  94. *     Jack Dongarra, Argonne National Lab.
  95. *     Jeremy Du Croz, Nag Central Office.
  96. *     Sven Hammarling, Nag Central Office.
  97. *     Richard Hanson, Sandia National Labs.
  98. *
  99. *
  100. *     .. Parameters ..
  101.       DOUBLE PRECISION   ZERO
  102.       PARAMETER        ( ZERO = 0.0D+0 )
  103. *     .. Local Scalars ..
  104.       DOUBLE PRECISION   TEMP1, TEMP2
  105.       INTEGER            I, INFO, IX, IY, J, JX, JY, KX, KY
  106. *     .. External Functions ..
  107.       LOGICAL            LSAME
  108.       EXTERNAL           LSAME
  109. *     .. External Subroutines ..
  110.       EXTERNAL           XERBLA
  111. *     .. Intrinsic Functions ..
  112.       INTRINSIC          MAX
  113. *     ..
  114. *     .. Executable Statements ..
  115. *
  116. *     Test the input parameters.
  117. *
  118.       INFO = 0
  119.       IF     ( .NOT.LSAME( UPLO, 'U' ).AND.
  120.      $         .NOT.LSAME( UPLO, 'L' )      )THEN
  121.          INFO = 1
  122.       ELSE IF( N.LT.0 )THEN
  123.          INFO = 2
  124.       ELSE IF( INCX.EQ.0 )THEN
  125.          INFO = 5
  126.       ELSE IF( INCY.EQ.0 )THEN
  127.          INFO = 7
  128.       ELSE IF( LDA.LT.MAX( 1, N ) )THEN
  129.          INFO = 9
  130.       END IF
  131.       IF( INFO.NE.0 )THEN
  132.          CALL XERBLA( 'DSYR2 ', INFO )
  133.          RETURN
  134.       END IF
  135. *
  136. *     Quick return if possible.
  137. *
  138.       IF( ( N.EQ.0 ).OR.( ALPHA.EQ.ZERO ) )
  139.      $   RETURN
  140. *
  141. *     Set up the start points in X and Y if the increments are not both
  142. *     unity.
  143. *
  144.       IF( ( INCX.NE.1 ).OR.( INCY.NE.1 ) )THEN
  145.          IF( INCX.GT.0 )THEN
  146.             KX = 1
  147.          ELSE
  148.             KX = 1 - ( N - 1 )*INCX
  149.          END IF
  150.          IF( INCY.GT.0 )THEN
  151.             KY = 1
  152.          ELSE
  153.             KY = 1 - ( N - 1 )*INCY
  154.          END IF
  155.          JX = KX
  156.          JY = KY
  157.       END IF
  158. *
  159. *     Start the operations. In this version the elements of A are
  160. *     accessed sequentially with one pass through the triangular part
  161. *     of A.
  162. *
  163.       IF( LSAME( UPLO, 'U' ) )THEN
  164. *
  165. *        Form  A  when A is stored in the upper triangle.
  166. *
  167.          IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN
  168.             DO 20, J = 1, N
  169.                IF( ( X( J ).NE.ZERO ).OR.( Y( J ).NE.ZERO ) )THEN
  170.                   TEMP1 = ALPHA*Y( J )
  171.                   TEMP2 = ALPHA*X( J )
  172.                   DO 10, I = 1, J
  173.                      A( I, J ) = A( I, J ) + X( I )*TEMP1 + Y( I )*TEMP2
  174.    10             CONTINUE
  175.                END IF
  176.    20       CONTINUE
  177.          ELSE
  178.             DO 40, J = 1, N
  179.                IF( ( X( JX ).NE.ZERO ).OR.( Y( JY ).NE.ZERO ) )THEN
  180.                   TEMP1 = ALPHA*Y( JY )
  181.                   TEMP2 = ALPHA*X( JX )
  182.                   IX    = KX
  183.                   IY    = KY
  184.                   DO 30, I = 1, J
  185.                      A( I, J ) = A( I, J ) + X( IX )*TEMP1
  186.      $                                     + Y( IY )*TEMP2
  187.                      IX        = IX        + INCX
  188.                      IY        = IY        + INCY
  189.    30             CONTINUE
  190.                END IF
  191.                JX = JX + INCX
  192.                JY = JY + INCY
  193.    40       CONTINUE
  194.          END IF
  195.       ELSE
  196. *
  197. *        Form  A  when A is stored in the lower triangle.
  198. *
  199.          IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN
  200.             DO 60, J = 1, N
  201.                IF( ( X( J ).NE.ZERO ).OR.( Y( J ).NE.ZERO ) )THEN
  202.                   TEMP1 = ALPHA*Y( J )
  203.                   TEMP2 = ALPHA*X( J )
  204.                   DO 50, I = J, N
  205.                      A( I, J ) = A( I, J ) + X( I )*TEMP1 + Y( I )*TEMP2
  206.    50             CONTINUE
  207.                END IF
  208.    60       CONTINUE
  209.          ELSE
  210.             DO 80, J = 1, N
  211.                IF( ( X( JX ).NE.ZERO ).OR.( Y( JY ).NE.ZERO ) )THEN
  212.                   TEMP1 = ALPHA*Y( JY )
  213.                   TEMP2 = ALPHA*X( JX )
  214.                   IX    = JX
  215.                   IY    = JY
  216.                   DO 70, I = J, N
  217.                      A( I, J ) = A( I, J ) + X( IX )*TEMP1
  218.      $                                     + Y( IY )*TEMP2
  219.                      IX        = IX        + INCX
  220.                      IY        = IY        + INCY
  221.    70             CONTINUE
  222.                END IF
  223.                JX = JX + INCX
  224.                JY = JY + INCY
  225.    80       CONTINUE
  226.          END IF
  227.       END IF
  228. *
  229.       RETURN
  230. *
  231. *     End of DSYR2 .
  232. *
  233.       END
  234.